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SUMMARY 

The convergence characteristics of various approximate factorizations for the 3D 
Euler and Navier-Stokes equations are examined using the von-Neumann stabil- 
ity analysis method. Three upwind-difference based factorizations and several 
central-difference based factorizations are considered for the Euler equations. In 
the upwind factorizations both the flux-vector splitting methods of Steger and 
Warming and van Leer are considered. Analysis of the Navier-Stokes equations 
is performed only on the Beam and Warming central-difference scheme. The 
range of CFL numbers over which each factorization is stable is presented for 
one-, two- and three-dimensional flow. Also presented for each factorization is 
the CFL number at which the maximum eigenvalue is minimized, for all Fourier 
components, as well as for the high frequency range only. The latter is useful 
for predicting the effectiveness of multigrid procedures with these schemes as 
smoothers. Further, local mode analysis is performed to test the suitability of 
using a uniform flow field in the stability analysis. Some inconsistencies in the 
results from previous analyses are resolved. 


INTRODUCTION 


Implicit numerical schemes are gaining increasing popularity since they allow 
large time steps for advancing the solution of Euler and Navier-Stokes equations 
to steady state. To reduce the computational cost that is usually involved, the 
implicit operator is often approximated by a number of smaller easily invertible 
factors. However, as observed by Thomas et al. [1], the approximately factored 
scheme has a stability restriction which is more severe in 3D, and also an optimal 
convergence time step that is not known a priori. Therefore, to avoid the long 
and costly approach of trial and error of obtaining an optimal CFL number, it is 
highly desirable to carry out a stability analysis for any numerical scheme. Some 
researchers have found that analyzing scalar equations such as the convection or 
the diffusion equation can provide insight into the stability requirements for Euler 
and Navier-Stokes equations. Beam and Warming [2] employed a combination 
of these scalar equations to approximate the restriction that will be placed on 
their ADI methods for compressible Navier-Stokes equations. Jameson and Yoon 
[3] and Caughey [4], among many others, used the scalar convection equation 
as a model problem for the Euler equations to investigate appropriate conditions 
for multigrid implementation. Rather than utilizing model equations, Jespersen 
and Pulliam [5] developed a technique where Fourier analysis is extended to the 
actual coupled equations of quasi-one-dimensional Euler equations. Jespersen [6] 
further extended this technique to 2D Euler equations in order to find the best 
conditions at which to implement multigrid for a transonic flow. Thomas et al. 
[1], von Lavante [7] and Anderson et al. [8] have also utilized a similar approach 
in the stability analysis of Euler equations for certain approximate factorizations 
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and relaxation schemes. 

In this paper, the stability analysis for both the Euler and Navier-Stokes equations 
is carried out for different approximate factorizations. For the Euler equations, 
three different upwind factorizations, the LU factorization and the Beam and 
Warming (ADI) factorization are considered while for the Navier-Stokes equa- 
tions, only the Beam and Warming (ADI) central scheme is analyzed. Also, the 
quasi-one-dimensional Euler equations investigated by Jespersen and Pulliam [5] 
is revisited in order to illuminate the actual influence on stability of using approx- 
imate Jacobians (to reduce computational costs), instead of the exact Jacobians 
in upwind factorizations. To ascertain the adequacy of using a uniform flow field 
in the stability analysis, a local mode analysis is further carried out using actual 
flow fields (transonic and subsonic) from a quasi-one-dimensional flow. 

THEORY AND ANALYSIS 


In order to extend the Fourier analysis to the coupled equations under consid- 
eration, a discrete analog of these equations is formulated based on different 
approximate factorizations in this section. The Euler equations are first analyzed 
using upwind and LU factorizations. The ADI factorization is formulated for the 
Navier-Stokes equations with the Euler equations as a degenerate case. 

Upwind Approximate Factorizations^ [or Euler Equations 

The conservation form of the 3— D Euler equations in Cartesian coordinates can 
be written as: 


dQ ,SE_ d^ dG = 

dt dx dy dz 

where Q is the solution vector and E, F and G are the 


( 1 ) 

conserved inviscid fluxes: 


3 



Q = [p,pu,pv,pw,pe] T 
E = [pu,pu 2 +p,puv,puw,(pe + p)u] 

F = [pv,puv,pv 2 + p,pvw,(pe+p)v] 

G = [pw,pwu,pwv,pw 2 +p, (pe + p)w] T 

If the Euler implicit scheme is used for time discretization, Eq. (1) can be written 
in the following form of the augmented Newton’s method [9]: 

[I + At(8 x A n + 8 y B n + 8 z C n )]AQ n = -A t{6 x E n + 8 y F n + 8 z G n ) (3) 

where the Jacobians A , B and C are and respectively. The 

expressions for A, B and C are given in Appendix A. 

Flux-vector splitting is employed for the upwind scheme where discretization of 
the flux derivatives is based on the physical propagation of the solutions of the 
Euler equations [10], Based on the direction of the characteristics at a grid point, 
A, B , C, E, F , etc., are split into their forward and backward contributions. 
Denoting the forward contribution with “+” and the backward contribution with 
and forward and backward difference operators with 8+ and 8~ respectively, 
we can rewrite Eq. (3) as: 

[1+ At(6;A+ + i+A') + Ai(«S-B+ + S*B-) + Ai{S;C* + 6tC~)iAQ = 

- A t{S;E + + S?E~ + S-F + + S+F~ + S;G* + StG~] 

The left hand side of the equation is usually approximated with first-order differ- 
ences, but the right hand side uses second-order differences to improve the overall 
accuracy of the converged solution. However, even with first-order difference ap- 
proximations of the implicit terms, the equation is computationally expensive to 
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solve. To reduce this cost, the implicit operator is factored into a sequence of eas- 
ily invertible terms. Following Anderson et al. [8] we will consider the following 
three factorizations: 

[i + Ai(6;A + + i+/4-)](i + a t(s;B + + s;b-)1 
[I + M{6~C+ + 6tC~))AQ n = -A tR n 
[i + A t(s;A + + s;B+ + s;c + )} 
p: + A t(8+A- + S+B- + 6+C~)]AQ n = -A tR n 

[i + At{s;A + + 6tA~ + s;c + )} 

(7) 

[I + A t(6~B + + 6~B~ -1- 6tC~)\AQ n = -A tR n 
Eq. (5), (6) and (7) shall be referred to as the spatial, eigenvalue and combination 

factorizations, respectively. 

There are different ways of obtaining the split fluxes expressed in the above 
equations but two popular methods viz: Steger and Warming flux-vector splitting 
[11], and van Leer flux-vector splitting [12], are considered in this work. 

In the Steger and Warming case, the fluxes are obtained from the following 
transformation: 


A + =X a D+X~\ A- = X a D~ a X~\ etc. (8) 


where D A and D A are diagonal matrices whose elements are the positive and 
negative eigenvalues of A, respectively, and the columns of X A are the eigenvec- 
tors of the Jacobian A. E + and E~ are obtained from E + = A + Q, E~ = A~Q 


etc. Eq. (8) gives approximate values for A + ,A etc. while exact values can 
be obtained from: 


8E+ _ dE~ 

~ ~dQ ’ ” dQ 


(9) 
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In order to resolve the singular nature of the Steger and Wanning flux-vector 
splitting at the sonic speed, a, van Leer proposed the following splitting in 
Cartesian coordinates: 


ri 


?± _ _^p( u + a ) 


E* = ± ! 


4 a 


((7 - l)u ± 2a) h 


( 10 ) 


w 


((7 - l)u ± 2a) 2 /2(7 2 - 1) + o ( y2 + ™ 2 ) 


With similar forms for F + , F~ ,G + ,G ~ , the Jacobians j 4 + ,.A _ etc. are obtained 
from Eq. (9). The analytical expressions for these can be obtained using a 
symbolic manipulator such as Mathematica. In these expressions, van Leer 
ensured continuous differentiability of the fluxes especially at the sonic transition 
[ 10 ]. 

LU Approximate Factorization (or Euler Equation^ 

This approach has become popular in recent times. It factors the implicit term of 
Eq. (3) into two components such that each component is strictly either a lower 
(L) or an upper (U) matrix as in the following equation: 

[I + MiS-Ai + S-Bi + £JCi)][I+ At{S}A 2 + S+B 2 + StC 2 )]AQ n = 

— At[8 x E 4- SyF 8 Z G) 

The Jacobian matrices are split to ensure diagonal dominance for each matrix 
inversion at each grid point. For our numerical computation we have adopted the 
flux-vector splitting devised by Jameson and Turkel [13]. 


* = <^I>, ^ = ^2. etc. 


( 12 ) 
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where r A > max(|Ayi|), etc. and are the eigenvalues of matrix A viz: 
u + a, u — a, u, u, u. 

The explicit terms are central differenced and it is necessary to damp the associated 
high frequency waves and/or to correct the odd-even decouplings. In this study, 
the following combination of second- and fourth-order explicit linear dissipations 
is employed. According to Caughey [4], and Yokota and Caughey [14], the former 
term is necessary for any spurious waves at the vicinity of shock while the latter 
ensures convergence to steady state. 

D e x = K2AtAx6 xx — K^AtAx 3 6 xxxx (13) 

Noting that 8 XX = ^-(<5+ -8~), the second-order term is split in a manner 
consistent with the differencing of the Jacobians and is implemented implicitly. 
Thus, with similar terms in the y- and z-directions, we write: 

[I + A t{8 x A l + S-Bi + Ci) + K 2 &t(6- + 8; + 8 7 )] 

[I + A t{8tA 2 + 8+B 2 + 8+C2) - K 2 At(8+ + 8+ + ^lAQ" = (14) 

-At{6 x E + 8 y F + 8 Z G ) - K 4 At(Ax 3 8 XIXX + A y 3 S yyyy + Az z 8 ZZIZ )lQ 

This factorization is similar to the eigenvalue factorization (see Eq. (6)) except that 
the explicit terms are centrally differenced rather than upwinded, thus, requiring 
the addition of dissipation. Also, the split fluxes of Jameson and Turkel which 
are less difficult to derive are used to achieve diagonal dominance in this case. 
ADI Factorization s for Euler_ and Na vier-S t okes Equations. 


The 3— D Navier-Stokes equations in Cartesian coordinates can be written as: 



where E, F and G are as defined earlier, and E v , F v , and G v are the viscous 
fluxes: 


E v — 


0 

|^(2u z - Vy- w z ) 
n(u y + V X ) 

p(u z + w x ) 

pv(u y + V x ) + pw(u z + W x ) + 
^pu(2u x — v y — w z ) + kT x 


(16) 


F v = 


0 

fi(Uy + V X ) 

\p{2v y -U x - w z ) 
fi(v z + Wy ) 

fiu(u y + V x ) + fiw(v z + Wy) + 

| [iv(2v y — u x — W z ) + kTy . 


(17) 


G v = 


0 

p(w x + u z ) 

li(V z + Wy) 

\p{2 W z - Vy - u x ) 
fiu(w x + Uz) + pv(v z + Wy) + 

\pw(2w z — Vy — u x ) + kT z _ 


(18) 


Where T = — A-tt and p is as defined in Appendix A. Also, Stokes hypothesis 
( A = -Ip.) has been assumed. With E v , F v and G v set to zero, we recover 
the Euler Eqs. (1). 

Following the approach of Beam and Warming [9], the viscous fluxes are split 
directionally. Also following the approach presented in Anderson et al. [14] 
for 2D Navier-Stokes equations, analysis yields the following ADI approximate 
factorization for the 3D Navier-Stokes equations, while assuming Euler implicit 


8 


time integration and constant fluid properties: 

[I + A t{6 x A - S XX R)][I + A t(S y B - S„S)]\l + A t{8 z C - 8 zz Y)]AQ n = 

— At[A8 x — R8 XX — R\8 yx ~ Rl^zx + Btiy ~ — S&yy ~ S 2 8 zy -{- 

C8 z -Y l 8 xz -Y 2 8 yz -Y8 zz )Q 

(19) 

The analytical expressions for the various Jacohians (from the viscous fluxes) that 
appear in this equation are shown in Appendix B. The right-hand-side resulted 
from linearization and from assuming the flux Jacobians locally constant. 

To damp the high frequency waves that will arise due to central differencing, 
second-order implicit ( D\ = — e,Af Ax<5 XI ) and fourth-order explicit 
(Dl = -e e AtAx z 8 xxxx ) artificial dissipations are added in the numerical exam- 
ples. Thus, with similar dissipations added in the y- and z-directions Eq. (19) 
becomes: 

[I + A t(8 x A - 6 XX R ) - e,-A<Axtfx*I][I + At{6 y B - 6 yy S) - e.AtAySyyl] 

[I + A t{6 z C - 8 ZZ Y) - £iAtAz8 zz I]AQ n 

— A t[A8 x — R8 XX — Ri8 yx — R 2 8 Z x + B8 y — S\8 xy — S8 yy — S 2 8 zy + C8 Z 

— Y\8 XZ — Y 2 8 yz — Y 8 ZZ + (seAx 3 8 XXX x + £eAy 8 yyyy + c e Az 8 ZZZZ )1]Q 


The corresponding factorization for the Euler equations is obtained by setting to 
zero the viscous flux Jacobians R,Ri,R 2 ,S,Si,S 2 ,Y,Y\,Y 2 . 

In the forgone analyses, different approximate factorizations that are widely used 
in practice have been formulated for the 3D Euler and Navier-Stokes equations. 
The convergence characteristics of each of these are examined using the von- 
Neumann type Fourier analysis methods. 
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von-Neumann Stability Analysis 

Each of Eqs. (5), (6), (7), (14) and (20) can be expressed as 

NAQ n = -L = —AtR n (21) 

von-Neumann stability analysis is used on this system of linear Eq. (21) by letting 
the step by step solution be characterized by 

Q n = A B [/oe 7 ‘^e 7 ^e 7 *^ (22) 

where I = y/— 1 , A is the amplification factor and <f> x , <j> y , <j> z represent the modes 
in the x-, y- and z-directions. Thus, Eq. (14) reduces to a complex generalized 
eigenvalue problem of the form [5]: 

I<v = ANv where K = N — L (23) 

The Fourier symbols jVand L are derived for each of the factorizations shown 
in Eq. (5), (6), (7), (14) and (20). For example, for the spatial factorization 
(represented by Eq (5)), employing a first-order differencing for the implicit 
operator and second-order differencing for the explicit operator, these two Fourier 
symbols are expressed as follows: 

N = jl + [(^ + _ ~ cos <?*) + (^ + + A~)I si n ^x] | 

jl + [(5+ - B~)( 1 - cos <p y ) + ( B + + B~)l sin <f> y ] J (24) 

jl+ ^[(C + - C")(l - cos + [C + + C~)I sin <f> z ] j 

L = — t r~\(A + - A - ) (3 + cos 2 <f> z - 4cos 4> x ) + (A + + A _ )(4sin<^ - sin 2 <j> x )l] 

2 Ax 1 

+ — [ (B + - B~) (3 + cos 2<f>y - 4 cos 4> y ) + ( B + + B~){ 4 sin <f> y - sin 2<f> y )l] 

2 Ay 

+ — f(C+ - C")(3 + cos2o, - 4 cos 6.) + (C + + C~)(4sine. - sin 2*,)/] 
2Az (25) 
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The Fourier symbols corresponding to the other approximate factorizations are 
documented in Demuren and Ibraheem [16]. 


SOLUTION PROCEDURE 


The convergence characteristics for solution algorithms based on each of the 
factorizations discussed are investigated by solving the generalized eigenvalue 
problem (23) over a fixed number of Fourier modes. 16 modes are selected, 
in the range 0 < <t> x ,4>yi<t>z < 2 tt, and over these modes the maximum eigen- 
value (Amax), the average eigenvalue (A avg ) and the smoothing factor (A^) are 
computed. The smoothing factor is computed to show the effectiveness of the 
selected scheme as a relaxation operator in a multigrid implementation. This 
is calculated from A^ = max(|A|) for the high frequency modes in the range 
f < <f>x, 4>y, <f>z < For the analyses, uniform flow is assumed with Moo = 0.8, 
zero yaw and angle of attack and 7 = 1.4. Further, the grid spacing is assumed 
to be uniform in all directions. The time step, A t is calculated from: 


At = 


CFL 




A* 


+ a V^ T + 




(26) 


As a further test case, quasi-one-dimensional Euler equations are solved with 
a similar formulation as the 3-D upwind spatial factorization, with uniform 
conditions of = 0 . 5 , zero yaw and angle of attack and p = 1 . 0 , chosen 
to enable comparison with Jespersen and Pulliam’s results [5]. In this case, the 
computed parameters are the maximum eigenvalue (A max ), the L2— norm of the 
eigenvalue (I 2 ) and the eigenvalue at d x = tt (A*-). 
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RESULTS AND DISCUSSIONS 

Computed values of the maximum eigenvalue (A max ), the average eigenvalue 
(A avg ) and the smoothing factor (A^) for the spatial, eigenvalue and combination 
factorizations based on the Steger and Warming flux-vector splitting are shown 
in Figs, (la), (lb) and (lc) respectively. Both the eigenvalue and the combina- 
tion factorizations are unconditionally stable for all CFL numbers. The spatial 
factorization is stable only for CFL numbers below 5. The maximum eigenvalue 
for each of the spatial, eigenvalue and combination factorizations is minimized at 
CFL numbers of 3, 8 and 7, respectively. Corresponding results obtained for 2-D 
case (not shown) indicate that the spatial and eigenvalue factorizations are uncon- 
ditionally stable and have lower (A max ) than the 3D case, for all CFL numbers. 
The corresponding minimum value of (A max ) are minimized at a CFL numbers 
of 8 and 10, respectively. The 1-D case is also stable for all CFL numbers with 
the maximum eigenvalue minimized at a CFL number of 11, for both spatial and 
eigenvalue factorizations (Table 1). 

Figs. (2a), (2b) and (2c) show the convergence characteristics of each of the 
factorizations based on the van Leer flux-vector splitting. These agree very 
well with that of Anderson et al. [8]. Except for the spatial factorization, 
all the schemes are unconditionally stable for all CFL numbers. The spatial 
factorization is stable only for CFL number below 14. The maximum eigenvalues 
for the spatial, eigenvalue and combination factorizations are minimized at CFL 
numbers of 7, 4 and 7 respectively. From the A^ curve, it appears that the 
spatial factorization with the Steger and Warming method has poorer smoothing 
properties comparison with the van Leer spatial factorization. Based on linear 
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analysis, there is also a smaller range of CFL numbers over which it is stable. The 
spatial factorization and the eigenvalue factorization of the 2-D case are found to 
be unconditionally stable with maximum eigenvalue minimized at CFL numbers 
of about 9 and 6, respectively. Results for the 1-D case are almost identical to 
those of the Steger and Warming analysis, with maximum eigenvalues minimized 
at CFL numbers of 11 and 19, respectively. 

In the computations presented thus far, approximate Jacobians derived from a 
time linearization of the Euler equations have been employed in the Steger and 
Warming method on both the implicit and explicit sides. The effect of using the 
exact Jacobians in the stability analysis was investigated with ID Euler equations 
using uniform conditions of Moo — 0.5 and p = 1.0. The results are compared 
in Figs. (3a) and (3b), respectively. In both cases, first-order differencing were 
used on the implicit side and second-order differencing on the explicit side, as in 
previous computations. From these figures, it can be observed that the results (as 
reflected by the variation of A max ,A a v g! A M with CFL) are similar. This shows 
that the use of approximate Jacobians does not place a restriction on the stability. 
This is at variance with the conclusion of Jespersen and Pulliam [5]. Restriction 
on the stability will result if the Jacobians are “mixed” such that approximate 
Jacobians are used on the implicit side and the exact Jacobians on the explicit 
side. In this case, Fig. (3c) shows that the stability is restricted to CFL numbers 
below 1. On the other hand, if the Jacobians are mixed in the reverse order 
i.e., with exact Jacobians on the implicit side and approximate Jacobians on the 
explicit side, the results (see Fig. (3d)) is not significantly affected. Further, from 
Figs. (4a), (4b), (4c) and (4d), where we have used second-order differencing on 
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both sides, similar conclusions can be drawn. 

All computations have been based on uniform how conditions. To ascertain the 
suitability of using such uniform flow field assumptions in the stability analysis, 
computations were carried out on two non-uniform flow fields with quasi-ID Euler 
equations using local mode analysis. These correspond to supersonic and transonic 
flows in a converging duct with steady-state solutions shown in Figs. (5a) and (5b), 
respectively. The von-Neumann method is applied at each point in the flow field 
thereby accounting for the variation in flow properties. The stability results for the 
supersonic case for both first-order and second-order differencing of the implicit 
side are shown in Figs. (5c) and (5d). Corresponding results for the transonic 
case are shown in Figs. (5e) and (5f). These results follow a similar trend as 
those obtained for ID Euler equations with uniform flow properties, except that 
instability is now predicted for lower CFL numbers. Boundary conditions were 
implemented explicitly and might have contributed to this instability. The use of 
local mode analysis here, is similar to the use of the total matrix method approach 
of Jespersen and Pulliam [5], except that, the former is easier to compute because 
it involves the solution of only a 3X3 eigenvalue problem. 

Figs. (6a), (6b) and (6c) show the convergence characteristics of the 3D Euler 
equations using the LU approximate factorization with central difference approxi- 
mations and various levels of second- and fourth-order artificial viscosities; « 2 and 
«4. Without' the addition of second-order dissipation i.e., /c 2 = 0, the coefficient 
/c 4 = 0.4 yields the optimal results (see Fig. (6a)). Appropriate combinations 
of K2 and *4 (especially, when *4 > ko) considerably reduce the amplification 
factor (see Fig. (6b) as compared with Fig. (6c)). The amplification factor is 
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minimized in each case at a CFL of about 5. Similar trends were observed in 
ID and 2D cases. 

In Figs. (7a), (7b), (7c), (8a), (8b) and (8c), the convergence characteristics for 
the full 3-D Naiver-Stokes equations using the Beam and Warming (ADI) central 
difference scheme as the baseline solution algorithm are shown for different 
Reynolds numbers and levels of artificial dissipation. For Reynolds number 
of 100 (Fig. 7a) and with no dissipation added, the scheme is stable for CFL 
number below 18. However, with artificial dissipation coefficients of e e = 0.5 
and £,• as 1.0 (Fig. 7b), the stability is restricted to a lower CFL number of 10, but 
with better smoothing properties. Optimal dissipation coefficients of s e = 1.0 and 
e,- = 2.0 (Fig. 7c), are found to improve the stability to a CFL of about 18 while 
maintaining good smoothing properties. The maximum eigenvalue is minimized 
at a CFL number of about 4 for this optimal dissipation. Both 1— D and 2— D 
cases are unconditionally stable for all levels of dissipation. For £ e = L0 and 
£,• = 2.0, their maximum eigenvalues are both minimized at about CFL numbers 
of 24 and 11, respectively. For Reynolds number of 10 6 , the results are similar 
to the cases with Reynolds number of 100, especially when dissipation is added. 
Hence, the stability results are not significantly affected by Reynolds number. 
Figs. (9a), (9b) and (9c) show the stability results for Euler equations with the 
Beam and Warming (ADI) central difference scheme. These results are identical 
to those obtained for the full Navier-Stokes equations at a Reynolds number of 
10 6 . Generally, the addition of dissipation reduces the amplification factor and 
the smoothing factor at lower CFL numbers. Optimal smoothing is usually at a 
CFL number close to 1. 
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The above results are summarised in Table I. In the Table, A ra stands for the 
minimum amplification factor, CFLm for the corresponding CFL number, CFL| 
the maximum CFL number for stability and CFL^ is the CFL number at which 
A n is minimized. 

CONCLUSIONS 

The stability of some approximate factorization schemes for the solution of the 
3D Euler equations and Navier-Stokes equations have been studied. For the 
Euler equations, the Steger and Warming, and van Leer flux-vector splittings 
were used with three different upwind factorizations namely: spatial, eigenvalue 
and combination factorizations. For both flux-vector splittings, the eigenvalue and 
combination factorizations are unconditionally stable, but the spatial factorization 
is only conditionally stable for CFL numbers below 5 for the Steger and Warming 
scheme, and 14 for the van Leer scheme. Moreover, the amplification factor 
(A m ax) minimized for the Steger and Warming scheme at CFL numbers of 3, 
7, and 8 respectively, and for the van Leer scheme at 7, 4, and 7, for spatial, 
eigenvalue and combination factorizations, respectively. Each of the approximate 
factorization methods has good smoothing properties for the van Leer flux-vector 
splitting, while for the Steger and Warming splitting, the smoothing factors are 
comparatively worse. Therefore, the van Leer splitting will be preferable for 
multigrid implementation. The Euler equations have also been analyzed for 
stability using the LU approximate factorization with central differences and 
various levels of artificial dissipation. It was found to be unconditionally stable 
in all dimensions with the maximum eigenvalue minimized at a CFL number of 
about 3. Contrary to the conclusion drawn by Jespersen and Pulliam [5] that the 
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use of approximate Jacobians places restriction on the stability, it is shown, after 
careful investigation, that if they are used on both the implicit and the explicit 
sides, the stability results are comparable to the case where the exact Jacobians 
are used. The von-Neumann analysis method was also employed in performing 
local mode analysis for actual (supersonic and transonic) flow fields of a quasi 
ID problem to show the suitability of using uniform flow field in the stability 
analysis. Stability results for the 3D Euler and Navier-Stokes equations solved 
with the Beam and Warming (ADI) central scheme with various levels of artificial 
dissipation (and at different Reynolds number tor the latter) have been presented. 
It was observed that the stability is not significantly affected by Reynolds numbers 
and that addition of dissipation reduces the amplification factor and the smoothing 
factor at lower CFL numbers. 
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APPENDIX A 


Inviscid Flux Jacobians 
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APPENDIX B 


Viscous Flux Jacobians 


R = 


dEv,i 

dQx 


a 

p 


s = 


dEj^y _ 

dQy 


Y = 


gg^ = 

dQ z 


R 1 = 


dEv, 

dQy 


y _ 




0 


0 


0 


0 

0 ' 





4 

3 


0 


0 

0 



— V 


0 


1 


0 

0 



—w 


0 


0 


1 

0 

*(9 2 - 

e ) 

1 40 2 

- ib ~ v - 

-u, 2 ) 

«(t- 

ib) 

V ( 1 - ib) 

w 1 

l 1 - ik) 

ib- 



0 


0 


0 


0 

0 ' 



—u 


1 


0 


0 

0 



~b 


0 


4 

3 


0 

0 



—w 


0 


0 


1 

0 


e ) 

- (b 2 - u 2 - 

w 2 ) 

u(l — 

P?) 


u>( 

1 - ib) 

ib . 



0 


0 


0 


0 

0 ' 



— u 


1 


0 


0 

0 



— V 


0 


1 


0 

0 



— iw 


0 


0 


4 

, 3 

0 

IO 

1 

• e ) 

A 0 0 

- (!“’ - U- 

-v 2 ) 

u(l — 

Pf) 

V ( 1 - pf) 

w( 

'4 7 \ 

k3 ~ P 7) 

0 

P? J 

' 0 

0 

0 0 O' 




' 0 0 

0 

0 O' 


0 

-f u 

0 

0 0 


dE v<z 


b o 

0 

-1 0 


— u 

1 

0 0 0 

,R 2 

~ dQ z 

— c 

p 

0 0 

0 

0 0 


0 

0 

0 0 0 



— U 1 

0 

0 0 



V 

-fti 0 0. 




b w w 

0 

-b o. 



dF ViX 

Sl 1 q7 ~ { 


0 


—7 


— U 
2 
3 
0 


0 0 0 

1 0 0 

0 0 0 

0 ° 0 0 0 

■|u u 0 0 


So = 


dF Vt2 

dQz 


0 0 
0 0 


0 
0 

b 

— V 

— | uv 0 w 


0 

0 


0 0 - 

0 1 0 


0 
0 
0 

0 0 

-b 0 


20 



y - dG *>* = t 

1 dQ x p 


0 

—w 

0 

f U 

— jUW 


0 0 0 0 

0 0 10 

0 0 0 0 

-§000 
-\w 0 u 0 


,^2 = 


9G V ,y 

dQy 


0 

0 

— w 
lv 


0 

0 

0 

0 


— juio 0 


0 O' 
0 0 

1 0 

-§ 0 0 
— §io v 0 


0 

0 

0 


where Pr = ^ 


21 


REFERENCES 


[1] J. L. Thomas, B. van Leer and R. W. Walters, Implicit Flux-Split Scheme For 
the Euler Equations, A1AA paper 85-1680, 1985. 

[2] R. M. Beam and R. F. Warming, An Implicit Scheme for the Compressible 
Navier-Stokes Equations, A1AA J., vol. 16. pp. 393—402, 1978. 

£3] A. Jameson and S. Yoon, Multigrid Solution ot the Euler Equations using 
Implicit Schemes, AIAA paper 85—0293, 1 985. 

[4] D. A. Caughey, Diagonal Implicit Multigrid Algorithm for the Euler Equations, 
AIAA J., vol. 26, pp. 841-850, 1988. 

[5] D. C. Jespersen and T. H. Pulliam, Flux Vector Splitting and Approximate 
Newton Methods, AIAA Sixth CFD Conference, 1983. 

[6] D. C. Jespersen, Design and Implementation of a Multigrid Code for the Euler 
Equations, Applied Mathematics and Computation, vol. 13, pp. 357-374, 1983. 

[7] E. von Lavante, D. Claes and W. K. Anderson, The Effects of various Implicit 
operators on a Flux Vector Splitting Method. AIAA paper 86—0273, 1986. 

[8] W. K. Anderson, J. L. Thomas and D. L. Whitfield, Three-Dimensional 
Multigrid Algorithms for the Flux-Split Euler Equations, NASA Technical Paper 
2829, 1988. 

[9] C. A. J. Fletcher, Computational Techniques For Fluid Dynamics 2, chap. 14, 
Springer- Verlag, New York, 1991. 

[10] C. Hirsch, Numerical Computation of Internal and External Flows 2, chap. 
20, John Wiley & Sons, New York, 1990. 


22 



[11] J. L. Steger and R. F. Warming, Flux Vector Splitting of the Inviscid Gas- 
dynamic Equations with Application to Finite-Difference Methods, J. of Comp. 
Physics, vol. 40, pp. 263-293, 1980. 

[12] B. van Leer, Flux-Vector Splitting For the Euler Equations, ICASE Report 
No. 82-30, 1982. 

[13] A. Jameson and E. Turkel, Implicit Schemes and LU Decompositions, 
Mathematics of Computaion, vol. 37, pp. 385-397, 1981. 

[14] J. W. Yokota and A. D. Caughey, LU Implicit Multigrid Algorithm for the 
Three-Dimensional Euler Equations, A1AA J., vol. 26, pp. 1061-1069, 1988. 

[15] D. A. Anderson, J. C. Tannehill and R. H. Pletcher, Computational Fluid 
Mechanics and Heat Transfer, chap. 9, McGraw Hill, New York, 1984. 

[16] A. O. Demuren and S. O. Ibraheem, Convergence Acceleration of the Proteus 
Computer Code With Multigrid Methods, Interim Report Prepared for the Internal 
Fluid Mechanics Division, NASA Lewis Research Center, 1992. 


23 



TABLE I: STABILITY ANALYSIS RESULTS FOR VARIOUS FACTORIZATIONS 
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Results for the ADI Euler Equations are identical. 
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(c) Combination Factorization 


Figs, (la)— (lc): Eigenvalue for Linear Stability 

Analysis of 3— D Euler Equations 

(Steger if Warming’s Upwind Factorization) 
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(c) Combination Factorization 

Figs. (2a) -(2c): Eigenvalues for Linear Stability 
Analysis of 3— D Euler Equations 
(van Leer’s Upwind Factorization) 
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(a) Approx. Jacobians (lhs, rhs) 
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(b) Exact Jacobians (lhs, rhs) 


Figs. (3a)-(3b): Eigenvalues for Linear Stability- 
Analysis of 1-D Euler Equations 
(First-order lhs; second— order rhs) 
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Figs. (3c)-(3d): Eigenvalues for Linear Stability- 
Analysis of 1-D Euler Equations 
(First-order lhs; second-order rhs) 
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(a) Approx. Jacobians (lhs, rhs) 
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(b) Exact Jacobians (lhs, rhs) 


Figs. (4a)— (4b): Eigenvalues for Linear Stability 
Analysis of 1-D Euler Equations 
(second— order both sides) 
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(d) Exact Jacobians lhs; Approx. Jacobians rhs 


Figs. (4c)— (4d): Eigenvalues for Linear Stability- 
Analysis of 1— D Euler Equations 
(second— order both sides) 
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(a) Supersonic case 
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(b) Transonic case 


Figs. (5a)— (5b): Steady solution of Quasi— ID Euler Equations 
(Pressure variation) 
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(c) First-order lhs; second— order rhs 
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(d) Second-order both sides 


Figs. (5c)-(5d): Eigenvalues for Linear Stability 
Analysis of 1— D Euler Equations 
(Local Mode Analysis, Supersonic case) 
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(e) First-order lhs; second— order rhs 
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(f) Second— order both sides 


Figs. (5e)-(5f): Eigenvalues for Linear Stability 
Analysis of 1— D Euler Equations 
(Local Mode Analysis, Transonic case) 







Figs. (6a)— (6c): Eigenvalues for Linear Stability 
Analysis of 3— D Euler Equations 
(LU Factorization) 
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(b) e,=0.5, e t =l 



Figs. (7a)- (7c): Eigenvalues for Linear Stability 
Analysis of 3-D Navier-Stokes Equations 
Beam and Warming’s ADI Factorization: Re=le2 
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Figs. (8a) -(8c): Eigenvalues for Linear Stability 
Analysis of 3— D Navier— Stokes Equations 
Beam and Warming’s ADI Factorization: Re=le6 
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Figs. (9a)-(9c): Eigenvalues for Linear Stability 
Analysis of 3— D Euler Equations 
Beam and Warming’s ADI Factorization 
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